#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import math

def square_root_period(n):
    v = []
    square_root = int(n**0.5)
    a, b, c = 1, 0, 1
    while True:
        r = (a * square_root + b) // c
        v.append(r)
        if r == 2 * square_root:
            return v

        b -= c * r
        m = math.gcd(math.gcd(a, c), math.gcd(b, c))
        a //= m
        b //= m
        c //= m

        x = a * c
        y = -b * c
        z = a * a * n - b * b
        a, b, c = x, y, z


def pell_equation(n):
    v = square_root_period(n)
    k = len(v) - 1
    for i in range(1, k + 1):
        v.append(v[i])

    p0, q0 = v[0], 1
    p1, q1 = v[1] * v[0] + 1, v[1]
    k = k if k % 2 == 0 else 2 * k
    for i in range(2, k):
        p2 = v[i] * p1 + p0
        q2 = v[i] * q1 + q0
        p0, q0, p1, q1 = p1, q1, p2, q2

    return p1, q1


if __name__ == "__main__":
    max = 0
    print("x^2 - D * y^2 = 1")
    for d in range(2, 1001):
        if int(d**0.5) == d**0.5:
            continue

        p, q = pell_equation(d)
        if p > max:
            max = p
            print(f"D = {d}, x = {p}, y = {q}")
